{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Access to ICESat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook is for providing codes for downloading and accessing to ICESat data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic concept of ICESat\n",
    "\n",
    "It does not provide water product (ATL13) unlike ICESat2, but it provides GLAH14 which shows land surface. According to literature reviews, they extracted the levels using water mask.\n",
    "\n",
    "I found 1) a py file for downloading ICESat directly.. and 2) a package I need to figure out how to use it.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Py file from NSIDC"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Package nsidc-subsetter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is a package to provide us to download ICESat 1 and 2, and ICEBridge: https://github.com/tsutterley/nsidc-subsetter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!git clone https://github.com/choms516/nsidc-subsetter.git"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This codes does not provide GLAH14... I am thinking of adding GLAH14 in this file if it is possible.\n",
    "\n",
    "First, just use its code to see whether it works well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cd ICESat_water_level/extraction/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cd nsidc-subsetter/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import codes.nsidc_subset_altimetry as ns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ff = './ICESat1'\n",
    "prod = 'GLAH12'\n",
    "vs = 34\n",
    "us_id = 'whaudtlr516@gmail.com'\n",
    "bound_box = [-50.33333,68.56667,-49.33333,69.56667]\n",
    "times = ['2009-01-01T00:00:00','2009-12-31T23:59:59']\n",
    "fmt = 'NetCDF4' "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ns.nsidc_subset_altimetry(ff, prod, vs, USER=us_id,\n",
    "    BBOX=bound_box, TIME=times, FORMAT=fmt, VERBOSE=True, UNZIP=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Failed today..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Open ICESat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "from pathlib import Path\n",
    "import h5py\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make a list of downloaded files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "##### load files\n",
    "## set the directory\n",
    "data_home = Path('/home/jovyan/ICESat_water_level/extraction/icesat/')\n",
    "## list them up and check them\n",
    "files= list(data_home.glob('*.H5'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are 89 files: which are within Tonle Sap Lake"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "89"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(list(files))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just open an file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "test = files[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<HDF5 file \"GLAH14_634_2103_002_0211_0_01_0001.H5\" (mode r)>\n"
     ]
    }
   ],
   "source": [
    "f = h5py.File(test, 'r')\n",
    "print(f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are 5 groups: ['ANCILLARY_DATA', 'BROWSE', 'Data_1HZ', 'Data_40HZ', 'METADATA']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ANCILLARY_DATA\n",
      "[]\n",
      " \n",
      "BROWSE\n",
      "['Image_00', 'Image_01', 'Image_02', 'Image_03']\n",
      " \n",
      "Data_1HZ\n",
      "['Time', 'Geolocation', 'Packet_data', 'Quality', 'Transmit_Energy', 'Reflectivity', 'Elevation_Flags', 'Atmosphere', 'DS_UTCTime_1']\n",
      " \n",
      "Data_40HZ\n",
      "['Time', 'Geolocation', 'Elevation_Surfaces', 'Elevation_Corrections', 'Elevation_Angles', 'Elevation_Offsets', 'Quality', 'Elevation_Flags', 'Transmit_Energy', 'Geophysical', 'Reflectivity', 'Waveform', 'Atmosphere', 'DS_DEMhiresArElv', 'DS_PeakNumber', 'DS_UTCTime_40']\n",
      " \n",
      "METADATA\n",
      "['COLLECTIONMETADATA', 'INVENTORYMETADATA', 'PROVENANCE']\n",
      " \n"
     ]
    }
   ],
   "source": [
    "for grs in list(f):\n",
    "    print(grs)\n",
    "    print(list(f[grs]))\n",
    "    print(' ')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([17.875, 17.923, 17.921, ..., 56.454, 58.095, 57.457])"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f['Data_40HZ']['Elevation_Surfaces']['d_elev'][:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "ename": "AttributeError",
     "evalue": "'slice' object has no attribute 'encode'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-21-d5b09b246fe6>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mf\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'METADATA'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'COLLECTIONMETADATA'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'Temporal'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32mh5py/_objects.pyx\u001b[0m in \u001b[0;36mh5py._objects.with_phil.wrapper\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32mh5py/_objects.pyx\u001b[0m in \u001b[0;36mh5py._objects.with_phil.wrapper\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32m/srv/conda/envs/notebook/lib/python3.7/site-packages/h5py/_hl/group.py\u001b[0m in \u001b[0;36m__getitem__\u001b[0;34m(self, name)\u001b[0m\n\u001b[1;32m    262\u001b[0m                 \u001b[0;32mraise\u001b[0m \u001b[0mValueError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Invalid HDF5 object reference\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    263\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 264\u001b[0;31m             \u001b[0moid\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mh5o\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mopen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mid\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_e\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mname\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlapl\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_lapl\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    265\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    266\u001b[0m         \u001b[0motype\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mh5i\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_type\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0moid\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/srv/conda/envs/notebook/lib/python3.7/site-packages/h5py/_hl/base.py\u001b[0m in \u001b[0;36m_e\u001b[0;34m(self, name, lcpl)\u001b[0m\n\u001b[1;32m    135\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    136\u001b[0m             \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 137\u001b[0;31m                 \u001b[0mname\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mname\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mencode\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'ascii'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    138\u001b[0m                 \u001b[0mcoding\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mh5t\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mCSET_ASCII\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    139\u001b[0m             \u001b[0;32mexcept\u001b[0m \u001b[0mUnicodeEncodeError\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mAttributeError\u001b[0m: 'slice' object has no attribute 'encode'"
     ]
    }
   ],
   "source": [
    "f['METADATA']['COLLECTIONMETADATA']['Temporal'][:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[]"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "list(f['METADATA']['COLLECTIONMETADATA']['ECSCollection'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I need to find information on GLAH14, specifically differences between 'Data_1HZ', 'Data_40HZ'. Also, I need to subset the files within TSL boundaries."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In my understanding 'Data_40HZ' is the original information. Let's include useful parameters, then I will subset it by ROI.\n",
    "\n",
    "Groups\n",
    "- Data_40HZ/Geolocation\n",
    "`d_lat`: latitude of shots\n",
    "`d_lon`: longitude of shots\n",
    "\n",
    "- Data_40HZ/Elevation_Surfaces\n",
    "`d_elev`: elevation without any corrections\n",
    "\n",
    "- Data_40HZ/Elevation_Corrections\n",
    "`d_satElevCorr`: saturation elevation correction (correction to elevation for sturated waveforms.. To apply it to the elevations it must be added to the elevation estimates\n",
    "\n",
    "- Data_40HZ/Quality\n",
    "`sat_corr_flg`: the saturation correction flag - the possible quality of the elevation data\n",
    "`i_satNdx`:understanding of concerns on data quality from saturation effects\n",
    "\n",
    "- Data_40HZ/Geophysical\n",
    "- `d_DEM_elv`: DEM elevation at the footprint location from the SRTM30\n",
    "\n",
    "- ANCILLARY_DATA\n",
    "- `internal_time_delay_date`: time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def glah14_to_df(filename):    \n",
    "    f = h5py.File(filename, 'r')\n",
    "    lat = f['Data_40HZ']['Geolocation']['d_lat'][:]\n",
    "    lon = f['Data_40HZ']['Geolocation']['d_lon'][:]\n",
    "    elev = f['Data_40HZ']['Elevation_Surfaces']['d_elev'][:]\n",
    "    sec = f['Data_40HZ']['Elevation_Corrections']['d_satElevCorr'][:]\n",
    "    scf = f['Data_40HZ']['Quality']['sat_corr_flg'][:]\n",
    "    satndx = f['Data_40HZ']['Quality']['i_satNdx'][:]\n",
    "    dem = f['Data_40HZ']['Geophysical']['d_DEM_elv'][:]\n",
    "    \n",
    "    glah14_df = pd.DataFrame({'Latitude':lat,'Longitude':lon,'Elevation':elev,\n",
    "                            's_El_Corr':sec, 's_Corr_f':scf,'in_sat':satndx,\n",
    "                            'DEM':dem})\n",
    "    return glah14_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "          Latitude   Longitude  Elevation  s_El_Corr  s_Corr_f  in_sat  \\\n",
      "0        75.376446  298.446866     17.875        0.0         0       0   \n",
      "1        75.377939  298.445066     17.923        0.0         0       0   \n",
      "2        75.379429  298.443276     17.921        0.0         0       0   \n",
      "3        75.380913  298.441503     17.842        0.0         0       0   \n",
      "4        75.382396  298.439736     17.891        0.0         0       0   \n",
      "...            ...         ...        ...        ...       ...     ...   \n",
      "1082195  37.859500  334.629687     57.231        0.0         0       0   \n",
      "1082196  37.861056  334.629417     58.130        0.0         0       0   \n",
      "1082197  37.862613  334.629147     56.454        0.0         0       0   \n",
      "1082198  37.864169  334.628876     58.095        0.0         0       0   \n",
      "1082199  37.865726  334.628605     57.457        0.0         0       0   \n",
      "\n",
      "                   DEM  \n",
      "0        1.797693e+308  \n",
      "1        1.797693e+308  \n",
      "2        1.797693e+308  \n",
      "3        1.797693e+308  \n",
      "4        1.797693e+308  \n",
      "...                ...  \n",
      "1082195  1.797693e+308  \n",
      "1082196  1.797693e+308  \n",
      "1082197  1.797693e+308  \n",
      "1082198  1.797693e+308  \n",
      "1082199  1.797693e+308  \n",
      "\n",
      "[1082200 rows x 7 columns]\n"
     ]
    }
   ],
   "source": [
    "test_df = glah14_to_df(test)\n",
    "print(test_df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Subsetting to the ROI."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "### bounds of Tonle Sap Lake\n",
    "sp_ex = [103.643, 104.667, 12.375, 13.287]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_df_subset = test_df.loc[(test_df['Longitude']>=sp_ex[0]) \n",
    "                  & (test_df['Longitude']<=sp_ex[1])\n",
    "                  & (test_df['Latitude']>=sp_ex[2])\n",
    "                  & (test_df['Latitude']<=sp_ex[3])]"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:notebook] *",
   "language": "python",
   "name": "conda-env-notebook-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
